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FOREWORD 


The  present  MBR  code,  written  in  the  Mathcad  (version  2001)  language,  was  developed  at  the 
Indian  Head  Division,  Naval  Surface  Warfare  Center  (NSWC).  It  represents  an  extensive  revision  and 
reformulation  of  an  earlier  code  that  was  written  in  Fortran  by  the  author  and  Mr.  Patrick  O’Neal  at 
NSWC  (White  Oak  Laboratory)  circa  1988.  Both  codes  were  based  on  a  theory  developed  by  the  author 
during  the  period  of  1982  through  1984,  which  showed  how  arbitrary  prior  distributions  for  ordered 
response  data  could  be  represented  by  a  mixture  of  ordered  Dirichlet  distributions  and  described  how 
rapid  and  exact  calculations  of  the  posterior  marginal  distributions  could  be  achieved  by  means  of 
recursive  relationships.  A  report  on  the  theory  has  been  recently  published  as  a  companion  to  this  report 
and  is  included  among  the  references.  Over  the  years,  the  MBR  program  has  been  used  by  the  Navy  in  a 
number  of  significant  applications  concerning  the  vulnerability  of  Naval  structures  to  explosions.  A 
digital  copy  of  the  code  can  be  obtained  by  contacting  Mr.  Hennessey  at  hennesseyts@ih .  navy.mil 
or  sending  a  request  to  the  Warheads  Branch. 


Timothy  S.  Hennessey 
Manager,  Warheads  Branch 


Approved  and  released  by: 


Amy  J.  O’Donnell 
Director,  Underwater  Division 


Released  by: 


iNW^5  

Marc  S.  Magdinec 
Head,  Weapons  Department 
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INTRODUCTION 


The  MBR  (Monotone  Bayesian  Regression)  computer  program  computes  posterior  marginal 
distributions  for  binomial  response  probabilities  associated  with  a  set  of  ordered  stimulus  levels.  Usually 
these  are  values  of  a  physics-based  measure  of  the  test  environment  severity.  We  shall  call  it  a  response 
initiation  index  function  or,  alternatively,  a  generalized  stress.  Use  is  made  of  posterior  marginal 
distribution  functions  for  Ramsey’s  M-variate  ordered  Dirichlet  joint  prior  (Ramsey,  1972),  which  were 
first  published  by  Disch  (1981).  Because  these  functions  involve  multiple  nested  summations,  they  are 
exceedingly  difficult  to  calculate  directly.  MBR  makes  use  of  recursive  expressions  developed  by  the 
author  (McDonald,  2003)  that  make  it  possible  to  achieve  rapid  and  exact  evaluations.  Moreover,  a  joint 
prior  consisting  of  a  single  M-variate  ordered  Dirichlet  distribution  is  not  sufficiently  general  to  apply  to 
most  problems  of  interest.  Consequently,  MBR  employs  a  related  joint  prior  that  is  a  mixture  of  M- 
variate  ordered  Dirichlet  distributions  to  obtain  a  class  capable  of  representing  quite  arbitrary  and 
general  forms.  The  mixed  prior  is  constructed  by  having  an  expert  provide  a  set  of  three  percentile 
curves  (including  the  median)  that  are  functions  of  the  initiation  index  as  shown  in  Figure  1.  Usually  the 
curves  are  chosen  as  the  10th,  50th,  and  90th  percentiles  of  the  prior  marginal  distributions.  The  code 
then  assigns  parameters  to  the  mixed  Dirichlet  prior  to  match  the  percentiles  and  calculates  the  posterior 
marginal  distributions  using  the  response  data  in  accordance  with  Bayes’s  law.  An  earlier  Fortran 
version  of  the  code  was  written  in  1988  by  the  author  and  Mr.  Patrick  O’Neill  at  the  Naval  Surface 
Warfare  Center  (NSWC)  (White  Oak,  MD).  The  current  Mathcad  code  was  written  at  the  Indian  Flead 
Division,  NSWC  and  completed  in  August  of  1999. 


Figure  1.  Prior  Marginal  Densities  and  Percentile  Curves 
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The  MBR  code  (version  5)  is  listed  in  the  following  section  and  contains  both  text  portions  and 
mathematical  coding.  The  code  has  an  interactive  fonnat  that  instructs  the  code  user  to  provide 
information  concerning  his  or  her  prior  percentile  curves,  to  input  the  test  data,  and  to  set  program 
constants  that  control  the  execution.  It  processes  the  data,  provides  a  graphical  display  of  the  posterior 
distribution  percentile  curves,  and  calculates  the  80%  coverage  interval  for  a  90%  probability  of  failure. 

The  present  code  is  written  in  the  Mathcad  (version  2001)  language.  Mathcad  was  chosen  because 
it  combines  a  user-friendly  textbook-like  style  with  very  powerful  mathematical  algorithms, 
programming,  and  graphical  capabilities.  Execution  proceeds  from  left  to  right  and  from  top  to  bottom. 

Previous  Mathcad  versions  of  the  code  were  two-dimensional  Mathcad  work  sheets  with  the 
interactive  narrative  and  results  running  down  the  left-hand  pages  of  the  worksheet  and  the  supporting 
equations  spread  out  to  the  right.  This  two-dimensional  structure  made  the  earlier  Mathcad  codes 
unsuitable  for  publication  as  reports.  In  an  effort  to  remedy  this  problem,  MBR  version  5  is  organized 
vertically  in  the  left-hand  pages  only  and  includes  four  areas  of  supporting  code  that  may  be  hidden 
from  view  to  enhance  program  readability  or  expanded  to  show  programming  details.  The  program  with 
the  detailed  coding  areas  hidden  is  best  for  users  whose  primary  goal  is  to  analyze  data.  The  areas  are 
expanded  in  this  report  to  show  all  elements  of  the  program.  The  hideable  areas  are  bounded  by 
horizontal  lines.  The  areas  may  be  hidden  or  unhidden  by  double  clicking  on  one  of  the  lines.  Generally, 
the  code  within  a  hideable  area  is  directly  relevant  to  the  section  that  follows.  The  function  of  the  coding 
within  each  area  is  indicated  by  a  numbered  AREA  title  written  in  8-point  font  above  or  below  the 
horizontal  line,  such  as 

U  AREA  1.  INPUT  &  REPRESENTATION  OF  THE  PRIOR  PERCENTILE  CURVES  (CODE) 

which  is  found  at  the  top  of  page  4.  The  texts  and  titles  within  the  collapsible  areas  are  presented  in 
italics  to  make  them  distinguishable  from  the  main  text.  Each  collapsible  area  can  also  locked  and 
password  protected. 
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MBR  COMPUTER  PROGRAM  (FULLY  EXPANDED) 

Introduction  to  MBR 

MBR  is  a  computer  program  that  is  used  to  express  and  update  (from  test  data)  the  probability  p 
of  some  arbitrarily  defined  binary  response  as  a  function  of  an  appropriately  chosen  quantity,  which  we 
symbolize  by  T  .  T  is  a  function  of  the  test  conditions  and  is  referred  to  as  a  response  intiation  index 
or  a  generalized  stress.  The  model  requires  that  p  depend  uniquely  upon  T  and  that  p  increase  (or  be 
nondecreasing)  with  T  .  No  other  assumption  concerning  their  functional  relationship  is  made.  The 
statistical  approach  is  Bayesian.  As  p  is  unknown  for  all  values  of  T  (excepting  possibly  at  zero  and 
infinity),  it  is  regarded  as  a  random  variable.  The  distribution  of  p  at  a  given  value  of  T  indicates  the 
uncertainty  of  the  response  under  the  conditions  implied  by  T  .  Prior  to  the  examination  of  test  data, 
MBR  requires  descriptions  of  the  distribution  of  p  as  a  function  of  T  over  the  range  of  T  values  of 
interest.  These  are  provided  in  the  form  of  three  distribution  percentile  curves.  MBR  then  constructs  the 
joint  prior  distribution  and  combines  it  with  the  binomial  test  data  via  Bayes  theorem  to  obtain  posterior 
marginal  distributions  and  updated  percentile  curves  that  indicate  how  the  initial  prior  uncertainties 
should  be  revised  in  light  of  the  test  results. 


Input  and  Representation  of  the  Prior  Percentile  Curves 

MBR  reconstructs  prior  marginal  distributions  from  three  probability  percentile  curves  —  the 
50th  percentile,  and  lower  and  upper  percentiles  chosen  by  the  user.  Enter  now  the  values  of  the  lower 
and  upper  percentiles: 


Lower  :=  10  Upper  :=  90 

To  input  the  curves  from  files,  skip  to  Tabular  Entry  of  Prior  in  Area  2  .  To  draw  the  curves 
graphically,  indicate  in  the  following  matrix  theY  values  where  each  curve  crosses  the  p  values 
indicated  by  the  column  headings.  BotLine  and  TopLine  percentile  values  are  arbitrarily  set.  Zero 
indicates  lower  threshold  values. 


BotLine  :=  .  1 

TopLine 

. 

.9 

^  "%ile  Curve  \ 

Zero 

BotLine 

.5 

TopLine 

Graphlnput  := 

Upper 

0 

1.9 

3 

4.5 

50 

0 

2.8 

4 

5.6 

V  Lower 

0 

3.8 

5 

6.6 
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Ij  AREA  1.  INPUT  &  REPRESENTATION  OF  THE  PRIOR  PERCENTILE  CURVES  (CODE) 


Code  for  Input  and  Representation  of  the  Prior  Percentile  Curves 


Special  Functions  Xcal  &  Peal  and  XBurr  &  PBurr 


Xcal(p  ,C ,/) 


x  0^  C\, o 

a  <-  Cp  i 

b  <-  Q, 2 

y  Q,3 

vvm  <—  Cj ,  4 

pr  <— - 

1-P 

1 

z  <— - 

—  Y 

1  +  ' 

1 

ex  <— - 

a  +  6-z 

ex  In  <—  ex-ln(  pr) 
x  <—  1037  if  exin  >  85 
xMo  if  exin  <  -85 

x<-xo  +  vvm  -pr  '  otherwise 


Pcal(p ,  pmi  ,C  ,i) 


a  <—  Q,0 

b  <-  Q,  l 


y  <-  Q,2 

i  -  p 
pr  <— - 

P 


pmr 


1  -  /;/«/ 
pmi 


ps^j 


a+- 


V  I  +  pmr-pr 


t+pr1 
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XBurr  ( ps ,  A  ,n)  := 


pr  B' 


pr  T  - 


1  -  PS  o 
ps  o 
l  -  ps2 


ps  2 

gl  <-  -ln{pr  B) 
gu  <-  — r) 
for  is  1 ..  f? 

*  0.  <-  ^i,o 


l 

Wl  <r 

wm 

wu  < 

Hi  <- 

rl  <- 


~  Ai,  1  _  ^i',0 

^  t ,  2  _  ^  ,  0 

“  3  ~  A,0 

Will 

wu 

wm 

wl 


yl  < - ln(rl) 

yu  < - ln(ru) 

y  <—  .9 

aa  < - 1 

bb  <-  -l 
count  <—  0 

while  [(aa  <  0)  +  (aa  +  bb  <  0)]  ^  0 
count  <—  count  +  1 
return  "count  =  20"  if  count  —  20 


y  +  .1 

l 


1  +  pr  B 
1 


y 


y 


bb  <- 

aa  <- 

at  <—  aa 
bj  <—  bb 
y vt  <-  y 
w/nv  ,•  <— 


1  +  pr  j 
gif  _gj_ 
yu  yl 
zu  -  zl 


gu  ,, 

- bb-zu 

yu 


augment  (x  q,  augment  {a  ,  augment  (b ,  augment  (y v  ,wmv  )))) 
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PBurr(ps  ,A  ,n)  •= 


wL<- 


w  u 


1-ps  0 
ps  0 
l-  ps2 


ps  2 

yL<~  - ln(w  L ) 

y  U  < - ln{w  u) 

for  is  1 ..  n 

1  _  4,  1 


zm  <— 


Ai,  1 
—In 


\Ai,0‘zm  ) 


h 


U  ■ 


—In 


'1-4,2^ 


Dif  < 


\  Ai,2'zm  J 

hjj  h  L 

yu  yL 

y  <-  .9 

aa  < - 1 

bb  <-  1 

C0i«2?  <—  0 

while  aa  <  bb -0.09984 
COt«2f  <—  COi«2?  +  1 
return  "count  =  20"  if  count  =  20 

y  <-  y  +  .1 
wTerm  <— - - - 


bb 


\  +  W  L 
Dif 


1 


—  wTerm 


y 


1  +  w  u 

h  1 

aa  <— - bb -wTerm 

yL 

—  aa 


bj  <—  bb 
Y Vi  <-  Y 
(a  b  yv ) 

augment{a,  augment  {b  ,yv)) 
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Tabular  Entry  of  Prior 

Code  to  create  Graphlnput  array  from  files  rather  than  by  the  user  should  be  inserted  here. 


Graph  Input  Processing  and  Display  of  Prior  Percentile  Curves 


Zero  =  "Zero" 

nplot  •=  30  j  :=  1  ■■nplot~  1 


pvj  :=  .5- 


/  •  \\ 

J-K 

1  -  COS 

V 

V  nplotj) 

T 

ps  :=  submatrix  (Graphlnput  ,0,0, 2,4) 

A  :=  submatrix  ( Graphlnput ,  0 , 3 , 1 , 4) 
XC  ;=  XBurr(ps,A, 3) 


PU  :=  pv  p5o  :=  pv  pl  :=  pv 

xu  ■=  Xcal^pu  ,XC,lj  *50  :=  Xcal^pso  ,XC ,2^  xl  :=  Xcal^pL  .,XC , 3) 
curveU  :=  augment  {x  u ,  p (j ) 
curve50  :=  augment  [x 50 ,  p 50) 
curveL  :=  augment  [xp,  py) 

ngj  :=  rows  (curveU )  -  1  //50  :=  rows(curve50)  -  1  ny  :=  rows  (cur  veil )  -  1 

ju  :=  1  ••  nu  j 50  :=  1  ••  nu  JL  :=  1  ••  *L 

Range  :=  ceil(max  (curveL)  +  l)  Range  =  12 
kk  :=  0..  1  xxq  :=  0  xx  1  :=  Range 


ft  AREA  1.  INPUT  &  REPRESENTATION  OF  THE  PRIOR  PERCENTILE  CURVES  (CODE) 
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Display  of  Input  Prior  Percentile  Curves 


PRIOR  PERCENTILE  CURVES 


RESPONSE  IN  IT  I  A  TION  INDEX 


Input  of  Observational  and  Nonobservational  Data 

Test  results  are  provided  to  MBR  through  a  four-column  matrix  called  “TestData.”  The  first 
column  is  the  test  number  (order  is  arbitrary),  the  second  column  is  the  value  of  the  test,  the  third 
column  is  the  number  of  tests  conducted  at  that  level,  and  the  fourth  column  is  the  number  of  responses 
observed.  Enter  information  into  the  TestData  matrix  now. 


r  "Test  Number" 

"Res. Init. Index" 

"#Tests" 

"#Responses 

1 

1 

1 

0 

2 

2 

1 

0 

3 

3 

1 

0 

4 

4 

1 

0 

5 

5 

1 

1 

6 

6 

1 

1 

TestData  := 

7 

7 

1 

0 

8 

8 

1 

0 

9 

9 

1 

1 

10 

10 

1 

1 

11 

11 

1 

1 

12 

12 

1 

1 

13 

13 

1 

1 

l  14 

14 

1 

1 
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Nonobservationai  Data 

In  addition  to  calculating  the  posterior  percentiles  at  the  values  of  T  listed  in  the  TestData  matrix, 
MBR  will  also  calculate  the  posterior  percentiles  at  an  arbitrary  number  of  nonobservationai  values 
evenly  distributed  across  the  response  initiation  index  range.  These  values  are  useful  both  for  predictions 
and  display  of  the  posterior  percentile  curves. 

Input  the  number  of  nonobservationai  values  desired:  nnonobs  :=  5 

The  Data  matrix,  augmented  by  nonobservationai  values  and  placed  in  order  of  ascending 
T  values,  is  shown  below.  (Note:  Data  of  any  duplicated  T  entries  are  combined.) 

JJ  AREA  2.  DATA  MATRIX  SUBMITTED  FOR  PROCESSING  (CODE) 


Code  for  Processing  Data  Matrix 

Range  for  Insertion  of  Nonobservationai  Data 


X05  ■=  Xcal(.05,XC ,2)  X.95  ■=  Xcal(.95,XC ,2) 

X.05  =  2.476  x.95  =  6.269 


ObsData  := 


return  0  if  nnonohs  -  0 
x.95~x05 

xint  ~ 

fhionobs  +  1 
Mat0  o  "nonobs" 

Maty,  1  X.05  +  xint 

Miato  2  ^ —  t) 

Mato ,  3  <—  0 

return  Mat  if  nnonohs  =  l 
for  k  e  1 ..  nnonobs  ~  1 
Matfc  o  "nonobs" 

Matk,  i  <-  Matk-\ ,  \  +  Xint 
Matk ;  2  <—  0 
Matk ,  3  <—  o 
Mat 
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Data 


return  TestData  if  ObsData  -  0 
Stack  <—  stack  {TestData ,  ObsData ) 
ndat  rows  {Stack)  -  1 

D  <—  submatrix  (  StocVc ,  l ,  nciat  •>  0 , 3) 

D  <—  csort{D,  l) 
nn  <—  0 

for  keo..ndat-l 
tf  Dnn ,  l  —  Dk ,  l 

Dnn  ,0  ^  Dk  ,0  f  Dnn,  0  =  "obs" 
Dnn ,  2  Dnn  ,2  +  ^,2 
Dnn ,  3  Dnn  ,  3  +  £>£ ,  3 
otherwise 


nn  <— 

77/7 

+  1 

^7171  ,  0 

<- 

A:,0 

^7171 ,  1 

<- 

Dk ,  1 

Dnn ,  2 

<- 

Dk,  2 

Dnn ,  3 

<- 

At,  3 

stack  [( "Test  Number"  "Res.Init.Index"  "#Tests"  "#Responses"  ),£)] 


(M  is  number  of  response  initiation  indices  in  Data  matrix.) 


M  :=  rows  {Data)  -  1  M  =  19  i  :=  1..M  T  :=  Data 


,<i> 


TOL  :=  10' 


10 


Given  Xcal{p  ,XC  ,k)  -  u  p  >  0  p  <  1  Xfun{p  ,XC  ,k  ,u)  :=  Find{p) 


p  :=  0.5 

Plo.  ■  =  Xfun (p,XC, 3 , T ,■) 
Pmidj  ■=  Xfun (p,XC, 2 ,  Y,-) 
Phi.  :=  Xfun  (p , XC ,  1 , T/) 


10 


IHTR  2323 


Graph  below  shows  p\0,  pmjd,  and  pjp  values  calculated  from  percentile  curves  at  the  test  and 
nonobservational  T  values. 


PRIOR  PERCENTILE  CURVES 


f, 1  AREA  2.  DATA  MATRIX  SUBMITTED  FOR  PROCESSING  (CODE) 

Data  Matrix  Submitted  for  Processing 


0 

1 

2 

3 

0 

"Test  Number" 

'Res.lnit.lndex" 

"#T  ests" 

"#Responses" 

1 

1 

1 

2 

0 

2 

2 

2 

1 

0 

3 

3 

3 

1 

0 

4 

"nonobs" 

3.108 

0 

0 

5 

"nonobs" 

3.74 

0 

0 

6 

4 

4 

1 

0 

7 

"nonobs" 

4.372 

0 

0 

8 

5 

5 

1 

1 

9 

"nonobs" 

5.004 

0 

0 

10 

"nonobs" 

5.637 

0 

0 

11 

6 

6 

1 

1 

12 

7 

7 

1 

0 

13 

8 

8 

1 

0 

14 

9 

9 

1 

1 

15 

10 

10 

1 

1 
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Data  Checking  and  Preparation  of  Dirichlet  Prior 
Prior  Distribution  Accuracy 

MBR  approximates  the  prior  marginal  distributions  as  mixtures  of  J  Beta  distributions  and  the 
joint  prior  distribution  as  a  mixture  of  J  ordered  Dirichlet  distributions.  The  number  of  terms  in  the 
mixtures,  J,  affects  the  accuracy  of  these  representations.  Usually,  J  set  in  the  range  of  20  to  50,  gives 
sufficient  accuracy  for  most  problems. 

Set  J  now:  J  :=20 

JJ  AREA  3.  DATA  CHECKING  AND  PREPARATION  OF  DIRICHLET  PRIOR 


Code  for  Data  Checking  and  Preparation  of  Dirichlet  Prior 

MBR  now  checks  that  ordering  constraints  imposed  on  the  prior  mixture  distributions  are 
satisfied. 

Locations  of  mixed  distribution  kernels  are  determined  by  equating  their  means  to  unevenly 
spaced  fractiles  of  the  reconstructed  prior  marginals.  The  spacing  is  determined  by  a  beta  distribution 
that  is  symmetric  about  p  =  0.5.  The  parameter  k  determines  the  degree  of  concentration  of  points 
toward  the  tails  of  the  distributions.  This  fitting  technique  performs  better  than  the  previous  method  of 
equating  the  modes  to  the  fractiles  of  equally  spaced  probabilities.  The  value  of  p  is  chosen  according 
to  the  number  of  kernels  J  so  that  twice  the  kernel’s  standard  deviation  equals  1/J  of  the  span. 


P curves 


(  Lower 

l,  100 


Upper 

100 


(  0.A 


P curves  ~ 


0.5 

VO.9J 


PA  :=  augment  ( piC) ,  augment  ( 


Pmid  i  Phi 


)) 


PC  PBurr(pcurves  ,PA  ,m) 


j:=  U.J 


*5 mult  ■—  1 


f  A 

ap  \^2jj  'amult  —  0.025 


P  := 


y  ®  mult  J 


p  =  397 
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Pord .  ■=  pbeta 
J 


\  J  +  1 


A 

,  1  +  K,  1  +  K 

) 


pstarjj  :=  Pcal[por(i. ,  pmid. ,  PC ,  /j 
:=  J  ai,  j  :=  (P  +  2)  -pstarjj 


b>i ,  j  •—  P  +  2  ai,j 


pstarjj 


P 


pstarjj  M+i  :=  1 


OrderingCheck  := 


for  j  e  1..J 
for  iel..M 

return  ( "Error  (i,j)  =  "  /  j  )  z/  \xstarj  i-\  >  [istarj  i 
"OK" 


Compare  Dirichlet  Marginal  (Mixed  Beta)  at  Y  h  with  Input  \iStar  Values 


J 


pMixedBeta{i , p ,a ,b ,©  ,j)  :=  ©  j'-pbeta^p ,aij> 


7=1 


dMixedBeta(i,p,a,b,0  ,J)  :=  © j'-dbeta(p,aij',bij>) 


j'  =  1 


fo  0]  Zt 

/zzze  :=  nk  :=  100  k  :=  o..zz&  /?/?£  :=  —  =  l 

VI  1/  nk 


J  =  20 


Set: 


h  :=  2  Fmb  :=  pMixedBeta ( h ,  pp/(  ,  a,  b,0  , j) 
h ,  k 


13 


IHTR  2323 


f|  AREA  3.  DATA  CHECKING  AND  PREPARATION  OF  DIRICHLET  PRIOR 


OrderingCheck  =  "OK" 

If  OrderingCheck  is  not  "OK",  try  reducing  the  value  of  J. 


Calculation  of  the  Posterior  Marginal  Distributions 

The  percentiles  of  the  posterior  marginal  distribution  functions  are  now  calculated  by  the  MBR 
and  are  plotted  below  alongside  the  prior  distribution  function  percentile  curves. 

Calculation  time  for  this  example  problem  takes  approximately  30  seconds  on  an  Intel  1200-MHz 
computer. 
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Ij  AREA  4.  CALCULATION  OF  POSTERIOR  MARGINAL  DISTRIBUTIONS  AND  DISPLAY  OF  POSTERIOR  PERCENTILE  CURVES 


Code  for  Calculation  of  Posterior  Marginal  Distributions  and  Display  of  Posterior 
Percentile  Curves 

Posterior  Parameters: 

M  =  19 

nfi  :=  Data i  2 
nti  :=  Data j  ,  i 

list  :=  nti  ~  nfi 

ntmax  '■=  max(nt)  ntmax  =  2 


mappo  :=  0  mappi  :=//[/=  1,1,  mappi- \  +  (/  -  2)  • ntmax  +  1  ] 
mapqo  :=  0  mapqt  :=//[/=  1,1,  mapqi- 1  +  (M  -  i  +  1)  -ntmax  +  1  ] 


nfle. 

l 


E  "h 


k  =  1 


ns 


'Sei 


M 

^  nsk 
k  =  i 


k  =  1 


M 

nfgej  ■=  ^  nfk 

k  =  i 
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CP-=  far  j'ei.J 


J.mapp  i 


break  if  M  <  2 
T  <-  1 
for  /e 


if  nf  iex  >  1 


t  <-  r- 


(°i.y  +  z~  0 
(fl2,y+  /- 1) 


for  k  e  0..  /e 

C  n  <—  J 

^  j,mapp2+k 

nkl  <—  nf  ie  +  k 


T -combin[ns\  ,k) 


T  <-  T- 


(a\j  + 

(a2,y  +  n^) 


break  if  M  <  3 
for  /e3..M 
T  <r-  1 

for  le  \  ..nfie 


if  nf  ie  >  1 

z— l 


T  <-  T 


{ai-l,j  +  i~  Q 

(a».y'+  0 


for  k  e  0..  /e 


S«W?  <—  0 


Al<r-(k  -  tlSj- 1  o) 
A2  ^k  ns  je  j 


^42  <—  ^  ns  ie  J 

for  r  e  max(Al) ..  min(A2 ) 

<—  sum  +  combiners i-\  ,k  -  /')  C 


Pj,mapp  ,_i+r 


C  n  . ,  <—  T -sum 

Pjjnapp  ,+k 

nkl  <—  nf  je  +  k 
i— l 


T  <-T 


(ai-ij+  nkl) 
{ai,j  +  nkl ) 
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Ca  :=  far  jel..J 


C  qj,mapqM  1 

break  if  M  <  2 
T<-  1 


for  l&\..ns , 


lf  nsgeM  >  1 


(%w'+ 


for  k  e  0..  nf , 


r 

qj,mapq^M_^+k 

nk  <—  ws  „„  +  k 
SeM 


T -combin(nfff,  k) 


T  <—  T  ■ 


( bM,j+nk ) 
(pM-fi+nk) 


break  if  M  <  3 
for  h  g  \  ..M  —  2 
i  <—  M  -  h  -  l 
T<r~  1 

for  l <=  l..ns  no 
Sei+\ 


if  ns  ge.+i  >  1 


r  <—  r  • 


(*i,/+7-0 


for  k  g  0..  «/, 


SM777  <—  0 


Al<-(k-nfi+i  o) 

"/gs+O 

for  s  g  max(Al) ..  min(A2) 
sum  <—  sum  +  combin(nf i+uk  -  s)-C  qjpnapq^+s 

C  qj,mapq  ,+k 


nk  <—  77s  „„  +  k 

se;+l 


r  <—  r  • 
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Posterior  Distribution  Function  Calculation 


UF(i,p)  := 


SUM  <-  0 
for  j  el../ 

n  <-  1 

c  <-  «,  ,y+  btj 

r<~nf  le. 

I 

if  r>  0 

„  aiJ 

n  < - - 


for  t  e  1 ..  r  —  1  if  r>  1 

(“ij+  0 


(c+  0 


if  s>  l 


n  <-  n- 

c  <-  c  +  r 
s  <—  ns  ge 

if  s  >  0 

bij 

n  ^  n — - 

c 

for  t  e  1 ..  s  -  1 

{pi  :+  t) 

n  <-  n  — - 

(c+  0 

T  1 

swm  <—  0 

^  <-  t>i,j  +  nf  ie , 


C<-  ^  +  S 

k  max  <-  =  1,0, ns 

k'  max  <“  */(*'  =  M  ,0,nf  J 

for  k<EO..kmax 

r  <-  1 
C'  <-  C  +  a- 
sum'  <—  0 
for  k'e0..k'max 
sum'  <—  sum'  ... 


T  <-  r  • 


+  (-0  C  y,mapqi+k'T'Pbeta  (P-1* 

B+k' 

C'  +  k' 


sum  <—  sum  +  (-1)  -C  D 

v  7  rj,mapp.+k 

A  +  k 


T  <-  r  ~ 


C+  k 


SUM 


k,B  +  k') 
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Note  that  the  normalizing  constant  FNorm  is  the  same  regardless  of  the  value  of  i  chosen  in 
UF(i,  l).  Disch  (1981)  regarded  this  as  a  check  of  code  accuracy. 

FNorm  :=  UF( l ,  l)  FNorm  =  1.773 x  io-4 


F(i,p) 


UF(i,p ) 
FNorm 


Distribution  Function  inverse  calculated  by  divide  and  conquer 


Finverse(ii ,  frac) 


win  <—  .0001 
PL<~  0 
PR  <-  1 
ii 

p< - 

M 

count  <—  0 
stop  <—  0 
while  stop  -  0 

count  <—  count  +  1 
return  p  if  count  >  100 
G  <—  F(ii,  p) 
if  G  <  ( 1  -  win)  ■ frac 
P  L  P 

P  L  +  P  R 

P  < - 

2 

if  G  >  ( 1  +  win )  • frac 
P  R^~  P 

PL  +  PR 

P  < - 

2 

stop  <—  1  otherwise 
P 
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>7,0  :=  Finverse(i,pcurves J 

>7,1  :=  Finverse(i,p  curves^ 

n,2  ■=  Finverse ( i , p curves] 

u '  :=  submatrix  (f,\,  M ,  0 ,  o) 
r'  \=  submatrix (r,  l, M ,0,2) 


vlo  \- 

v50  := 
vup 


k  :=  0..50 


ujc  :=  —  -max  ( u ') 

/  :=  0..M  -  l 


0 

1 

2 

0 

0 

0 

0 

1 

0 

0 

3.71431049547697-10  -6 

2 

3.19376244748894-10  -10 

7.95063219572368-10  -4 

0.012669613486842 

3 

4.14276123046875-10  -3 

0.026810495476974 

0.129060444078947 

4 

6.10753109580592-10  -3 

0.034796463815789 

0.155530427631579 

5 

0.03620750025699 

0.128045333059211 

0.368215460526316 

6 

0.064978348581414 

0.196443256578947 

0.471808182565789 

7 

0.133983812834087 

0.325336657072368 

0.613898026315789 

8 

0.334189967105263 

0.574022795024671 

0.793778268914474 

9 

0.3359105963456 

0.575709292763158 

0.794793379934211 

10 

0.582815872995477 

0.769980982730263 

0.89464689555921 

11 

0.703163548519737 

0.842336554276316 

0.927477384868421 

12 

0.894301565069901 

0.942192479183799 

0.971880461040296 

13 

0.959650541606702 

0.976931120219984 

0.988926937705592 

14 

0.982188375372636 

0.991033453690378 

0.996671174701891 

15 

0.990755583110609 

0.996504131116365 

0.999245091488487 

16 

0.99473089920847 

0.998739945261102 

0.999904532181589 

17 

0.996898450349507 

0.99964061536287 

0.999996260592812 

18 

0.998219540244655 

0.999938262136359 

0.999999982059786 

19 

0.99906325340271 

0.999995857477188 

0.999999999997101 
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pIoj(  :=  inter p{  vlo  ,u' ,r'^  ,u^) 
p50fc  :=  inter p{v 50  ,u',r'^  ,u/( 

(  ( 2 ) 

pupk  '■=  inter p\vup  ,u' ,r'  ,u/( 

w'o  +  u'm-  1 

uu  := - 

2 


up  :=  root 


( interp{vup , u' , r'®  ,  uu)  -  .9 ,  uu)  up 


up  :=  root 


{interp{vlo  ,u' ,r'^  ,uu )  -  .9 ,uu) 


UR 


5.688 

7.053 


ft  AREA  4.  CALCULATION  OF  POSTERIOR  MARGINAL  DISTRIBUTIONS  AND  DISPLAY  OF  POSTERIOR  PERCENTILE  CURVES 

Posterior  Percentile  Curves 


POSTERIOR  &  PRIOR  PERCENTILE  CURVES 
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80  Percent  Coverage  Interval 

The  80%  coverage  interval  (ur,ur)  for  the  response  index  associated  with  a  0.9  probability  of 
response  is 


UL  =  5.68  ur  =  7.05 


Notes — Current  Code  is  Version  5 


MBR3:  Version  3  differs  from  Version  2  in  the  following  respects:  Version  3  calculates  the 
posterior  density  function  along  with  the  distribution  function  and  computes  the  distribution  function 
percentiles  by  a  Newton-Raphson  scheme.  (Version  2  used  a  divide-and-conquer  technique.)  Plots  of  the 
posterior  density  function  and  distribution  function  are  available.  (Look  below  the  FNorm  calculation.) 

MBR4:  Version  4  goes  back  to  the  divide-and-conquer  scheme  for  calculating  the  inverse 
posterior  distribution  function.  Although  faster,  the  Newton-Raphson  algorithm  was  found  to  be  less 
robust.  Posterior  density  function  coding  used  in  Version  3  remains. 

MBR5:  Version  5  is  the  code  of  this  document.  It  differs  from  Version  4  only  in  the  use  of  a 
vertical  fonnat  and  hideable,  lockable  areas. 
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CONCLUDING  COMMENTS 


MBR  appears  to  be  functioning  properly.  Future  improvements  contemplated  include  development 
of  a  more  efficient  algorithm  for  calculating  the  inverse  posterior  distribution  as  currently  performed  by 
Finverse.  Possible  choices  include  a  hybrid  Newton-Raphson  and  divide-and-conquer  scheme. 

The  PBurr  and  XBurr  functions  calculate  the  distribution  function  and  inverse,  respectively,  of  a 
modified  form  of  the  Burr  distribution.  The  Burr  distribution  can  be  found  discussed  in  Kendall  and 
Stuart  (1960).  Representations  were  needed  in  MBR  of  a  distribution  function  and  its  inverse  whose 
parameters  were  the  median  and  two  percentiles.  These  are  used  for  graphing  purposes  and,  most 
importantly,  to  construct  modified  Burr  representations  of  the  prior  marginals  from  the  percentile 
curves,  which  are  then  used  to  assign  parameters  to  the  mixed  beta  forms  of  the  prior  marginals  and  the 
mixed  Dirichlet  joint  prior.  Difficulties  may  be  encountered  with  the  modified  Burr  distributions  when 
the  prior  percentile  curves  are  unusually  or  improperly  positioned  relative  to  each  other.  Hence,  these 
too  may  be  the  subject  of  future  improvements. 
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